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Abstract 

Quantum systems with a non-conserved probability can be described by means of non-Hermitian 
Hamiltonians and non-unitary dynamics. In this paper, the case in which the degrees of freedom can 
be partitioned in two subsets with light and heavy masses is treated. A classical limit over the heavy 
coordinates is taken in order to embed the non-unitary dynamics of the subsystem in a classical 
environment. Such a classical environment, in turn, acts as an additional source of dissipation 
(or noise), beyond that represented by the non-unitary evolution. The non-Hermitian dynamics 
of a Heisenberg two-spin chain, with the spins independently coupled to harmonic oscillators, is 
considered in order to illustrate the formalism. 
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I. INTRODUCTION 


Historically, the development of non-Hermitian quantum mechanics started with the 
study of metastable states and tunneling by means of Hamiltonians with complex energy 


eigenvalues 


2|]. Later a connection to PT-symmetry was found [3| and the concept 


of pse„do-herm.t.city was also estabUshed fl Q. Nowadays, this field of teseateh is con- 
stantly growing. A hrst general book on the topic has appeared [6|; applications of non- 
Hermitian quantum mechanics involve the study of scatte ring by complex potentials and 
quantum transport jz !?!, description of metastable states 18-23|, optical waveguides Q- 
261, multi-photon ionization 27-29], and nano-photonic and plasmonic waveguides j^. The 


theoretical investigations are also undergoing rapid developments: non-Hermitian quantum 
mechanics has been investigated within a relativistic framework 


31| and it has been adopted 


32 -£|. 


Moreover, it 


by various researchers as a means to describe open quantum systems 
seems that a few theoretical studies have been dedicated to the statistical mechanics and 
dynamics of systems with non-Hermitian Hamiltonians ^ 


3 - 


In the present work, the interest is focused on the development of a formalism to embed 
consistently the quantum dynamics of systems with a non-conserved probability in a classical 
environment, which is explicitly taken into account (he., it is not averaged-over) in the dy¬ 
namics and which , in turn, acts as a source of disorder. Types of noise beyond those arising 
from Gaussian white noise ^ can then be treated. From a more general perspective, one 


goal of this work is to develop a numerical formalism (which is complementary to that based 
on master equations [5^) for studying the dissipative dynamics of, for example, quantum 


plasmonic metamaterials 


56 


57l | or processes of interest in quantum thermodynamics 


In order to obtain such a formalism, a composite system with heavy and light degrees of 


freedom is considered. First a classical limi 
using the partial Wigner representation 59 


over the heavy degrees of freedom is performed 


67| . In such a way, a general quantum-classical 


approximation of non-Hermitian quantum mechanics is obtained. Finally, the limiting case 
in which the non-Hermitian part of the evolution does not affect the classical-like degrees of 
freedom (represented in Wigner phase space |68|, l69|) is considered. 

This paper is structured as follows. In Sec. |TT] the quantum non-Hermitian equation of 
motion for the density matrix is taken as a starting point and the classical limit over the 
heavy mass coordinates is performed by using the partial Wigner transform and the linear 
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expansion in the sqnare root of the ratio between light and heavy masses. From this general 
case, one can easily derive the eqnation of motion valid for a decay operator depending only 
on the quantnm degrees of freedom of the snbsystem. In Sec. IIIII the Hermitian part of the 
total Hamiltonian of the system is considered; an adiabatic Hamiltonian is extracted from 
this and its eigenstates are nsed to represent the qnantnm-classical non-Hermitian eqnation. 
In Sec. HV] piecewise-deterministic algorithms (using the adiabatic basis) are presented ex¬ 
plicitly in the case when the decay operator depends only on the quantum coordinates of the 
subsystem. The numerical approach is illustrated in Sec. |V]by studying (in the case of a weak 
coupling to the environment) the dynamics of a chain of two spins, coupled separately to an 
independent harmonic oscillator. Two different decay operators are explicitly treated. The 
evolution of the trace of the reduced density matrix of the spin chain and of relevant matrix 
elements is monitored. The results show that non-Hermitian quantum-classical dynamics 
(and the numerical algorithms developed in this work) can model the loss of probability and 
the damping expected in open quantum system. Conclusions are hnally given in Sec. I VII 


II. NON-HERMITIAN QUANTUM MECHANICS IN A CLASSICAL BATH 


Consider a composite quantum system with quantum coordinates {f,p,R,P) = {x,X). 
A multidimensional notation will be adopted in the following so that, for example, R stands 
for (i?i, i? 2 , Rn), where N is the total number of degrees of freedom in the conhgurational 
space of the subsystem represented by the operators X. It is also assumed that the dynamics 
of the composite system is dehned by the non-Hermitian Hamiltonian 


H = H-it , 


( 1 ) 


where H = [R + 'W)/2 and f = i{R — 'w)/‘^ Hermitian operators. In particular, in 
order to develop the formalism, one can assume that 


H 


p2 

2M 



+ V{r,R) . 


( 2 ) 


while T is left unspecihed. In order to have an idea of the type of physical systems to 
which the formalism arising from Eq. ([1]) can be applied, one can consider, for example, a 
quantum resonance of the following type: a double well with a decay operator providing 
the metastability of the bound state. The non-unitary dynamics would allow the particle to 
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escape from the bound state while eventual quantum tunneling would allow it to go from 
one well to another. Such internal transitions between the wells can be enhanced by the 
thermal disorder provided by a classical environment coupled to the well. 

In the following, it will also be assumed that M, the mass associated to the sub-system 
with coordinates X, is much bigger than m, the mass associated to the subsystem with 
coordinates x, i.e., M » m. This leads to the dehnition of the small parameter /i = 
The non-normalized density matrix of the composite system with non- 
Hermitian Hamiltonian "H obeys the equation of motion ^ 

1 




- hi 




( 3 ) 


where and are the commutator and anticommutator, respectively. 

In order to obtain the quantum-classical limit of Eq. ([3]) a rigorous procedure, which is 
based on the partial Wigner representation of the dynamics and the linear expansion in 
/i, can be followed. Such a procedure was used for Hermitian Hamiltonians in Ref. 641. 
Accordingly, one can introduce the partial Wigner transform of H over only the coordinates 
of the heavy degrees of freedom: 


^ - zimm+ zm . ( 4 ) 

As a results VLY^{X,t) is an operator in terms of the quantum x variables and a function 
in terms of the X variables (which are still quantum - they have only been represented in 
the Wigner quantum phase space). Analogously, an arbitrary quantum operator x of fh® 
composite system is partially transformed in Wigner phase space as 


Xw(X) = j - Z/2\x\R + Zjl) . (5) 

Moreover, the partial Wigner transform of a product of arbitrary operators x ^-nd ^ is given 
by 

(xi) {X) ^ , ( 6 ) 

\ / w 

where V = {{d/dR), (d/dP)) is the phase space gradient operator and 

0 1 

B= (7) 

-1 0 
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is the symplectic matrix, so that is basically the Poisson bracket operator, which will 

also be denoted with the symbol 

Upon taking the partial Wigner transform of Eq. ([3]), one obtains 
^nw(X,f)= - 

( 8 ) 


Equation ([8]) is still fully quantum in nature. Now, in order to take the quantum classical 
limit one can follow the procedure of Ref. 6^, which was inspired by the theory of Brownian 


motion given in 


73|, and introduce scaled coordinates. Arbitrary units can be dehned 


m 


the following way: one can introduce cq as unit of energy, to = fi/eo as unit of time, and 
Ao = h /as unit of length. One can also define po = and Pq = (Afeo)^'^^ as 

the unit of the light and heavy momenta, respectively. As in 6^, adimensional coordinates, 
making the momenta of the same order of magnitude, can be introduced: 


r = 


Ao 


p' = ^ 
^ PO 


JD/ _ R 

“ Ao ’ 

p/ _ _P^ 

^ Po ■ 


( 9 ) 


Accordingly, one has H^{x,X) = eoH{^{x, X) and fw(5^,-A) = eofw(^)-A). In the scaled 
coordinates given in Eq. ([9]), Eq. (E]) becomes 






( 10 ) 


where t' = t/to- Now one can take advantage of the smallness of p and expand the expo- 
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nential operators retaining only the linear order terms. One obtains 






fW(A"),siW(A",0 


+ E {x'XXx'}) VMX'X) 

jk 

- ^B,t{v'fi'Xx'X})x'xi„(x') 

- |Bjt(v'fW(A'))v'A(A',i') 


( 11 ) 


where the snm over the indices j, k runs over all phase space dimensions. Equation flTT|) 
gives the quantum-classical approximation to Eqs. (jS]) and ([8]) in adimensional coordinates. 
Transforming back to fully dimensional variables, one hnally obtains 


m 


,t) 


h L 

+ 2 t)| 


1 
h L 


rw(x),Ow(x,f) 

t), 

' ^'fw(A),tiw(A,«)| + |nw(A,«),fw(A) 


( 12 ) 


Equation (IT^ provides the rigorous quantum-classical approximation for the non-Hermitian 
dynamics of composite systems. It is a valid approximation when the degrees of freedom of 
the system have two different De Broglie wavelengths, one short and one long. It is worth 
noting that, formally, Eq. flT^ could have been obtained directly from Eq. ([8]) by taking an 
expansion to linear order in the limit h —>■ 0, without going through the transformation to 
scaled coordinates given in Eq. ([9]). However, the /i expansion seems much more rigorous 
since h ^ 0 would imply a cancelation of quantum effects also on the x coordinates. Never¬ 
theless, the practical agreement of the two limiting procedures indicates that the linear form 
in Eq. flT^ can perhaps simply assumed as an ansatz for quantum-classical dynamics, as 
suggested in Ref. [^. It is also very important to remark that Eq. flT^ has a form that does 
not depend on any particular basis. Moreover, Eq. flT^ is exact when the decay operator 
rw(-^) is linear in X, Hxn{X) is at most quadratic in the X coordinates and is linearly 
coupled through the Xs with the quantum subsystem. 

While Eq. flT^ dehnes non-Hermitian quantum dynamics in a classical bath in the case 
of a general decay operator rw(W), there is one interesting limiting situation that can be 
considered. It concerns the case in which the decay operator does not depend on the bath 
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coordinates. When this happens the partial Wigner transform leaves f invariant, so that 
Eq. flT^ reduces to 


t) 


2 


r L J + 


(13) 


Equation flT^ shows that, in this case, the effects on the subsystem dynamics arise from the 
anticommutator of E and r2w(-^) alone. Equation flT^ might describe a situation in which 
two types of effects are present: the non-Hermitian dynamics, with its probability leakage or 
pumping, of a quantum subsystem embedded in a bath of classical degrees of freedom, whose 
influence is expressed through the Poisson bracket terms on the right hand side. If also the 
density matrix does not depend on the bath coordinates X, VL-^{X,t) —)■ one obtains 

the purely quantum case, given by Eq. ([3]). Such a result constitutes a self-consistency check 
for the formalism. 


III. REPRESENTATION IN THE ADIABATIC BASIS 


The partial Wigner transform of the Hamiltonian in Eq. ([2]) can be rewritten as 

ffw(A') = ^ + im{R ). (14) 

The adiabatic basis of H^{X) is defined by the eigenvalue problem 

hY4\a] R) = Ea{R)\a] R) . (15) 
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Such a basis can be used to represent Eq. ([12]). Upon defining the quantities 


= {a;R\h^^{X,t)\a';R) , 
Ea{R) - Ea'{R) 


UJ rvrv' 


JpOi _ 

-^w — 


h 

dE^{R) 


dR ’ 

,r ^ L A , , p-b A 

M dR^ W'* dP ’ 

d 

f P 

^afi I ■ di^fS j h^ajsdafi ; 

'Taa',P/3' = ' 5P 

+ A w* A + A* A 


^^OLo! “1“ ^-^aaO ’ 

'^I^OLo! ^^olol' 'Taa',j3l3' : 


(16) 

( 17 ) 

(18) 

(19) 

( 20 ) 

( 21 ) 


apy • 

(23) 

(24) 


and using the results of section III in Ref. [6^, providing the representation of the first, 
third and fourth terms in the right hand side of Eq. (IT^ . one obtains for Eq. (IT^ : 


dt 


',/3/3'f^w i'w,^w{X,t) ^\a';R) 


013 ' 

^{a;R\ (|fwAw(X,f)} + |A(X,f),fw}) \a'-R) . (25) 


To proceed one can take advantage of the identities 


-ia7 


(“i-Rl^b;-R) = ^ - E “ E rw<*.7. 


{7;R| 


dR 

t) I , 


dR 
d 


dR 


o';-R> = ^fii"'-E<7f*w'-E‘*-'f*w 


(26) 

(27) 










The derivation is lengthy but straightforward. Its details are provided in App. The 
representation of Eq. ffT3|) into the adiabatic basis is 



- 


'(Ta^ W 


a 



(28) 


Equation (128|1 provides the desired result: it is the representation of the quantum-classical 
non-Hermitian dynamics, given by Eq. flT^ . in the adiabatic basis of ifw(-A). As it stands, 
it looks quite formidable and not amenable of being numerically integrated in a simple way. 

When the decay operator depends only on the quantum coordinates of the subsystem, 
one has to consider the representation of Eq. (IT^ in the adiabatic basis. It can be verihed 
that Eq. fl25D reduces to 



(29) 


Equation fl2^ is much simpler than the general Eq. fl28|l : it can be used to model dissipative 
effects on quantum subsystems arising both from sources and sinks of probability and dis¬ 
order because of the coupling to a classical environment. Equation fl29|) can be integrated 
by means of piecewise-deterministic algorithms. These will be sketched in Sec. IIVI 
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IV. PIECEWISE-DETERMINISTIC ALGORITHMS 


Equation (|28|) . which is the adiabatic representation of the abstract Eq. (na, can be 
numerically integrated by means of piecewise-deterministic algorithms. However, given the 
complexity of Eq. (125]) . only the limiting case given by Eq. (12 9 p will be explicitly considered 
here. 

Consider Erst Eq. (I29p . In such a case, it is convenient to decompose the representation 
of the decay operator in the adiabatic basis in terms of a diagonal, r^“, and an off-diagonal 
part, 

ra/3 ^ ^ . (30) 


Equation (l2^ becomes 


m 


p/3' ^ PP' 

1 


-R E (rf !JW'(A', t). 


PP' 


At this point, it is useful to dehne 




PQ;Q _I_ pc 

H “T i ^ 


^a. a. 
d “T -L ^ 


h 


TL.m' = \ (ro“'’W + rf-'&j) 


Equation (l3Tp becomes 

d 


a «w' = - E + TL.m) sjy (A.t). 

PP' 


where one has dehned 


(31) 


(32) 

(33) 


(34) 


^^OLo! 3B' “ 1 “ '^olol' “ 1 “ ^0.(3^a' • 


(35) 


A piecewise-deterministic algorithm for the integration of Eq. (13^ can be found by using the 


sequential short-time propagation (SSTP) scheme [66|, A trajectory can be seen as the 
concatenation of small hnite time steps At. Accordingly, for a single step, the propagator 
associated with Eq. (IMl) can be written as 


,-iAt(£(T')+r+rr) 


—iAtC 


Ct) 


olol' , 1313 ' 


iPoL^^a.' p' ^t'Taa'13' ,^I3'^ 


(36) 
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The propagator decomposition in Eq. fl36|) can be used as the basis for a SSTP algorithm for 


integrating Eq. flSTll. ' 


using either basic 


66 


^he actions of Tap^a'iS' and TSaSn>fi> must be sampled probabilistically 


75| or more advanced schemes 76l-l78l| for efficient convergence. The 


momentum-jump approximation 




79] can be adopted in the expression of Taa' ,13/3'■ It 


should be noted that the damping (or enhancing) frequency is considered in the action 


of iC, 


(7) 


V. NON-HERMITIAN SPIN CHAIN IN HARMONIC BATHS 


Consider a subsystem described by the Hermitian Hamiltonian 




u:, 


3z(y\ 


Z ^ Z 5 


(37) 


which represents a chain of two coupled quantum spins, s^, k = 1,2. The constants ji, with 
i = x,y, z, determine the spin coupling strength. The operators are the given by the 
Pauli matrices for spin kg = 1,2. The excited and ground state of the spins are denoted by 
|e(sfc) > i^(sfc) y, — 1^2), respectively. As in Ref. 80(], the subsystem basis is defined 


by the following vectors: |1 >= |2 >= >, |3 >= >, 

|4 >= >. The bath is composed by two harmonic oscillators and has the following 

partially Wigner-transformed Hamiltonian: 




H-—-K, 


k=l 


2M 


(sk) 


(38) 


Equation fl38l) provides the Wigner-transformed Hamiltonian of two independent harmonic 
oscillators with mass M and frequency u. Oscillator 1 is coupled to spin 1 while oscillator 
2 is coupled to spin 2. The coupling Hamiltonian (in the partial Wigner representation) is 


H< 


SB,W — 


-Y^cR 


(Sfc) 


a. 


(Sfc) 


(39) 


k=l 


The total partially Wigner-transformed Hermitian Hamiltonian of the system is (of course) 
given by Hy/{X) = Hs + Hsb,w + Hb,w- Since the total bath is harmonic and the coupling 
with the spin chain is bilinear, the linear approximation of the partially Wigner represented 
dynamics is exact. This means that the classical-like representation of the bath in Wigner 
space is, in fact, fully quantum in nature. 


11 


















In order to illustrate the numerical implementation of the formalism, two decay operators 
are considered: 


f = 71/ , 


(40) 

(41) 


The symbol I denotes the identity operator in the Hilbert space of the spin chain while 7 ^, 
j = 1,2, are constants. The operators are chosen so that the difference between the effect of 
a uniform probability sink on all states of the spin-chain, represented by operator T^^^ and 
the depletion of just the one state (when both spins are excited), represented by operator 
can be observed. The dynamics of the density matrix is determined by substituting 
j = 1 ^ 2 , into Eq. (IT^ . The adiabatic basis representation of Eqs. (IT^ has been given 
in Sec. HU Non-adiabatic corrections to the dynamics can be disregarded upon assuming 
a weak coupling to the environment; This means that the transition operators Taa',/ 3 i 3 ' and 
Xya',i 3 i 3 ' are neglected in the calculations here discussed. If At is the numerical 

integration step, a single-step SSTP propagator in the adiabatic approximation is written 


as 




1 pAt 


(42) 


where 70 , 0 ' is dehned in Eq. fl32H . The right hand side of Eq. flT2D can be derived by means 


64| . In the calculations reported either or 


of the Dyson identity, as explained in Ref. 

T*^^^ have been used to obtain 'yaa'-, depending on the case. 

In order to perform the numerical study, the density matrices of the subsystem VL^{t) and 
of the oscillators f 2 B,w(A, f) are considered uncorrelated at the inital time 


Dw(A, h) — Ds(h) Db,w(A, t\) , 


(43) 


where 


VL 


B,W 


(Vi,) = n 


tanh(/5a;/2) 


ks=l 


TT 


exp 


otanh(/3a;/2) ' 

-/- ) 


bJ 


(44) 


with = l/k-QT inverse thermodynamics temperature {k-Q denotes the Boltzmann constant) 
and H-q^n^X) is dehned in Eq. (|38ll . The initial condition for the reduced density matrix of 
the spin chain has been chosen as Dg(ti) = |<I) >< <I)| with |<I> >= when using T*^^) 

and as Clsik) = IT >< T|, with |T >= (l/\/ 2 ) >) when using 
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FIG. 1: Adiabatic time evolution of the trace of the reduced density matrix, Tr[IIs], for j3 = 
0.1, jx = jy = —1, jz = 0.5, c = 0.24. The numerical time step of integration is At = 0.01. 
Adimensional parameters are used. Results for non-Hermitian dynamics with decay operator F^^) = 
7 ii and initial reduced density matrix Cls{to) = |4> >< <I>|, with |<I> >= >. The upper 

curve shows the results for 71 = 0 (Hermitian dynamics). Then, from top to bottom, the curves 
for 71 = 0.1,0.5,1 are displayed. All curves are drawn with statistical error bars. 

Upon choosing cu as a reference frequency, one can introduce a unit of energy, hu, a 
dimensionless time, t —)■ ut, and an inverse thermodynamical temperature, [3 —)■ Tiufi. One 
can also introduce the adimensional coordinates R {h/MojY^^R and P (huMy^'^P. 
Correspondingly, one can use the following adimensional constants ji —)■ ji/fioj {i = x, y, z), 
c —)■ c/{hMuy, ■jk —t Tfc/hcu {k = 1,2). The values adopted in the calculations have been 
(3 = 0.1, ix = iy = —1, iz = 0 . 5 , and c = 0.24. The parameters specifying the decay 
operators have been set, in different calculations, to 71 = 0.1, 0.5,1 and 72 = 0.001, 0.01, 0.1. 
Phase space averages with negligible statistical errors have been calculated using 5 x 10^ 
points. 

Figured] displays the adiabatic time evolution of the trace of the reduced density matrix 
of the spin chain, Tr)!!^], when the decay operator is T^^^ = 711 . and initial reduced density 
matrix fis(O) = |*F >< <F|, with |$ >= >. The upper curve shows the results 

for 7 i = 0 (Hermitian dynamics). Then, from top to bottom, the curves for 71 = 0.1, 0.5,1 
are displayed. All curves are drawn with statistical error bars (which are already negligible 
by using just 50000 phase space points). As expected the “loss of probability” (given by 
the fact that non-Hermitian dynamics represents in an effective way the effect of additional 
states, which do not appear in the Hamiltonian, whose occupation can grow at the expense 
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t 

FIG. 2 : Adiabatic time evolution of matrix element in the subsystem basis, for /3 = 0.1, ja; = 
jy = —1, jz = 0.5, c = 0.24. The numerical time step of integration is At = 0.01. Adimensional 
parameters are used. Results for non-Hermitian dynamics with decay operator F^^ = 7 ii and 
initial reduced density matrix Ag(to) = |4> >< 4>|, with |<I> >= >. The upper curve 

displays the results for 71 = 0 (Hermitian Dynamics). Then, from top to bottom, the curves for 
71 = 0.1,0.5,1 are displayed. Starting from the top curve, a constant shift of 1.5 in the negative 
y direction has been applied for visualization purposes. All curves are drawn with statistical error 
bars. 

of the occupation of the explicitly described states; such is the case when escaping from 
a well toward inhnity or decaying from a metastable state) increases upon increasing 71 . 
Figure [2] displays the damped time evolution of the reduced matrix element of the spin 
chain. The results confirm that the numerical algorithm is stable and that, at least for 
the model considered, it can be used to simulate the decay of quantum states in a classical 
environment. 

Figure [3] displays the adiabatic time evolution of the trace of the reduced density 
matrix of the spin chain, Tr[f 25 ], when the decay operator is F^^^ = 72 |e^®i\ 
g(«i)^g(« 2 )| initial reduced density matrix Gg(0) = |T >< T|, with IT >= 

(l/\/2) >). The upper curve shows the results for 72 = 0.001. 

Then, from top to bottom, the curves for 72 = 0.01 and 72 = 0.1 are shown. Figure 0 ] 
displays the corresponding damping phenomenon for the reduced matrix element Gg^ of the 
spin chain. In this case, the trace decays because of the depletion of the state >, 

as it can be verihed by monitoring the time evolution of the diagonal elements of the density 
matrix in the basis of the spin system. 
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FIG. 3; Adiabatic time evolution of the trace of the reduced density matrix, Tr[flg], for /3 = 
0.1, jx = jy = —1, jz = 0.5, c = 0.24. The numerical time step of integration is At = 0.01. 
Adimensional parameters are used. Results for non-Hermitian dynamics with decay operator F^^) = 
72 |e(^i), g(si)^g(s 2 )| initial reduced density matrix Cls{to) = I'h >< Tl, with |T >= 

(l/\/2) >). The curves, from top to bottom, show the results for 72 = 

0.001,0.01,0.1. All curves are drawn with statistical error bars. 

VI. CONCLUSIONS 

In this work, a formalism to embed non-Hermitian quantum dynamics in a classical bath 
has been provided. In order to achieve this, a quantum-classical approximation for the 
non-Hermitian equations of motion of composite systems (with degrees of freedom having 
light and heavy masses, m and M, respectively) has been hrst been considered, using a 
partial Wigner representation. Then, the limiting case when the non-Hermitian part of the 
evolution does not involve the classical-like degrees of freedom has been taken into account. 
The classical bath embedding the quantum system with non-conserved probability can be 
used as a noise source that is more general than those of Gaussian type, e.g, colored noise. 
The bath can be also used to describe thermal disorder. 

When the adiabatic part of the Hermitian Hamiltonian is considered, its eigenstates 
(dehning the adiabatic basis) can be used to represent the non-Hermitian quantum-classical 
equation of motion. Once the equations of motion are represented in this adiabatic basis, 
algorithms can be developed using a sequential short-time propagation scheme. For the sake 
of illustrating the formalism, a Heisenberg chain with two spins, each weakly coupled to a 
separate harmonic oscillator has been studied. Two different decay operators have been 
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FIG. 4: Adiabatic time evolution of matrix element in the subsystem basis, for f5 = 0.1, 
jx = jy = jz = 0.5, c = 0.24. The numerical time step of integration is At = 0.01. 
Adimensional parameters are used. Results for non-Hermitian dynamics with decay operator 
f(2) = e(si)^e(« 2 )| initial reduced density matrix lls(to) = >< 'h|, with 

|T >= (l/\/2) >). The curves, from top to bottom, show the results 

for 72 = 0.001,0.01,0.1. All curves are drawn with statistical error bars. 

explicitly considered showing that the algorithms lead to a stable and efficient numerical 
approach. 

Future applications will be devoted to the modeling of nano-scale solid state devices in 
dissipative environments. 

Acknowledgements 

This work is based upon research supported by the National Research Foundation of 
South Africa. 

Appendix A: Representation of the equation of motion 

Starting from Eq. fl25l) the derivation can proceed by considering the term 

(a; R\ [fw, fiw(X, t)] + \a'; R) = (rw W + K "'• 

/ 3 / 3 ' ^ 

(Al) 
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Using Eq. flAip . the equation of motion in fl2^ becomes 


' = - E ^ E (^W t) 


/3/3' 


/3/3' 


— -(a; R\ |fw, ^)| 1“^ R) — -^1 “! ^), fw 1“ W] R) 


One can then consider the third term in the right hand side of Eq. flA2p : 


(A2) 


(a;fl|{fw.f!w(.’f,i)||a';-R> = ^ | (a; fl|^|7; fl) 




dP 


dR 


(A3) 


Using the identities in Eqs. (12611271) . Eq. flA3ll becomes 

r '1 / (9 

{a-R\ {fw,Ow(X,t)} 


HP’ 


_Ka __ 




E 


/3/3' L o- ^ 

d 

dP 


+ r^(io-/3^a'/3'77TA 


a^'(x,t) 


+ 


/3/3' \ a- 


w 


dP 


dP 


(A4) 
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Using Eq. flA4j) . Eq. flA2D becomes 



cr 


— -(a;/2| ^), fw| -R) • 


(AS) 


In order to complete the derivation, one must consider the last term in the right hand side 
of Eq. (IASI) : 
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dR 


^ elf'll / 


dP 


-b;R\- 


dR 


dP 

a'] R) 


(A6) 


7 

To further simplify Eq. flASp . one has to consider again (just with different indices) the 
identities given in Eqs. (12611271) . Using Eqs. (1261127p . Eq. (IA6P provides the hnal representation 
of the equation of motion in the adiabatic basis. 
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